Set Up

library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.3.2
## Warning: package 'tidyr' was built under R version 4.3.2
library(here)
library(phyloseq)
library(vegan)
## Warning: package 'lattice' was built under R version 4.3.2
theme_set(theme_bw())
filt_fun <- function (x, min_reads = 100, min_samples = 5) {
    (sum(x) > min_reads) & (sum(x > 0) > min_samples)
}
ps <- readRDS(here('data','following_study','ps_rarefied.rds')) %>% 
    filter_taxa(filt_fun, prune = TRUE) 
sample_data(ps)[is.na(sample_data(ps))] <- 'NA'
sam <- data.frame(sample_data(ps))
max.core <- parallel::detectCores()

Alpha Diversity

alpha.diversity <- data.frame(sample_data(ps),estimate_richness(ps),read_depth = rowSums(otu_table(ps)))

Household Effect

ggplot(alpha.diversity,aes(x = Household, y = Shannon, group = Household)) +
    geom_point() + geom_line()

anova(lm(Shannon~Household, data = alpha.diversity))
## Analysis of Variance Table
## 
## Response: Shannon
##           Df Sum Sq Mean Sq F value  Pr(>F)  
## Household 48 8.4359 0.17575   1.566 0.06064 .
## Residuals 49 5.4992 0.11223                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Epileptic Effect

ggplot(alpha.diversity,aes(x = Epileptic.or.Control, y = Shannon)) +
    geom_boxplot() + geom_jitter(height = 0, width = 0.25)

Paired t test

epileptic <- alpha.diversity %>% filter(Epileptic.or.Control == 'Epileptic')
control <- alpha.diversity %>% filter(Epileptic.or.Control == 'Control')
t.test(epileptic$Shannon, control$Shannon, paired = TRUE)
## 
##  Paired t-test
## 
## data:  epileptic$Shannon and control$Shannon
## t = -0.22662, df = 48, p-value = 0.8217
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.1529073  0.1219302
## sample estimates:
## mean difference 
##     -0.01548853

ANOVA with Household effect added

anova(lm(Shannon~Household + Epileptic.or.Control, data = alpha.diversity))
## Analysis of Variance Table
## 
## Response: Shannon
##                      Df Sum Sq  Mean Sq F value  Pr(>F)  
## Household            48 8.4359 0.175747  1.5357 0.07042 .
## Epileptic.or.Control  1 0.0059 0.005877  0.0514 0.82168  
## Residuals            48 5.4933 0.114444                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Breed Effect

ggplot(alpha.diversity) +
    geom_point(aes(x = Breed.Group..1., y = Shannon, colour = Breed.Group..1.)) +
    facet_wrap(~Epileptic.or.Control) + 
    theme(axis.text.x = element_blank(), axis.ticks.x.bottom = element_blank())

# here breed group with na is removed 
anova(lm(Shannon~Household + Breed.Group..1., data = alpha.diversity %>% filter(Breed.Group..1. != 'NA')))
## Analysis of Variance Table
## 
## Response: Shannon
##                 Df Sum Sq Mean Sq F value  Pr(>F)  
## Household       45 7.2456 0.16101  1.6626 0.05393 .
## Breed.Group..1.  4 1.3804 0.34511  3.5637 0.01431 *
## Residuals       39 3.7768 0.09684                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Drug Effect

alpha.diversity.epi <- alpha.diversity %>% filter(Epileptic.or.Control == 'Epileptic')
ggplot(alpha.diversity.epi, aes(x = Pheno.Y.N, y = Shannon)) +
    geom_boxplot() + geom_jitter(height = 0, width = 0.25)

t.test(Shannon~Pheno.Y.N, data = alpha.diversity.epi)
## 
##  Welch Two Sample t-test
## 
## data:  Shannon by Pheno.Y.N
## t = 0.84674, df = 33.509, p-value = 0.4031
## alternative hypothesis: true difference in means between group No and group Yes is not equal to 0
## 95 percent confidence interval:
##  -0.1450834  0.3521435
## sample estimates:
##  mean in group No mean in group Yes 
##          2.730402          2.626872

Sex Effect

ggplot(alpha.diversity, aes(x = Sex, y = Shannon)) +
    geom_boxplot() + geom_jitter(height = 0, width = 0.25)

anova(lm(Shannon~Household + Sex, data = alpha.diversity))
## Analysis of Variance Table
## 
## Response: Shannon
##           Df Sum Sq Mean Sq F value  Pr(>F)  
## Household 48 8.4359 0.17575  1.6253 0.04789 *
## Sex        1 0.3088 0.30882  2.8560 0.09752 .
## Residuals 48 5.1903 0.10813                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
t.test(Shannon~Sex, data = alpha.diversity)
## 
##  Welch Two Sample t-test
## 
## data:  Shannon by Sex
## t = -1.315, df = 87.559, p-value = 0.1919
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -0.25315537  0.05154705
## sample estimates:
## mean in group F mean in group M 
##        2.658970        2.759775

Age Effect

ggplot(alpha.diversity, aes(x = as.numeric(Age..months.), y = Shannon)) +
    geom_point()
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

lm(Shannon~as.numeric(Age..months.), data = alpha.diversity) %>% summary()
## Warning in eval(predvars, data, env): NAs introduced by coercion
## 
## Call:
## lm(formula = Shannon ~ as.numeric(Age..months.), data = alpha.diversity)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.92260 -0.28131 -0.01761  0.27021  0.78472 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              2.6317958  0.0896466  29.357   <2e-16 ***
## as.numeric(Age..months.) 0.0007684  0.0010859   0.708    0.481    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3742 on 94 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.005298,   Adjusted R-squared:  -0.005284 
## F-statistic: 0.5007 on 1 and 94 DF,  p-value: 0.4809

Beta Diversity

ord <- ps %>% ordinate(method = "MDS", distance = 'unifrac')
plot_ordination(ps, ord)

The PCoA plot shows a very separate pattern of four clusters. Let’s take a look at the phlogenetic tree of our data.

plot_tree(ps, "treeonly")

Some of the taxa have extremely long branch. Are they influential?

# Warning: Not all nodes are descendants in this tree
if (require(adephylo)) {dist.root <- adephylo::distRoot(phy_tree(ps))}
## Loading required package: adephylo
## Loading required package: ade4
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'RNeXML'
hist(dist.root); dist.root %>% sort(decreasing = TRUE) %>% head(20)

##    ASV601    ASV438    ASV563    ASV287    ASV484    ASV133    ASV562    ASV437 
## 20.819315  5.082379  4.338564  4.267155  4.245708  4.159508  4.154919  4.121791 
##    ASV450    ASV286    ASV478    ASV357    ASV446    ASV454    ASV423    ASV380 
##  4.088335  4.077083  4.074518  4.060563  4.045067  4.041199  4.026320  1.631610 
##    ASV213    ASV644    ASV268     ASV41 
##  1.617836  1.569113  1.417223  1.411706
detach('package:adephylo', character.only = TRUE, unload = TRUE, force = TRUE)

Let’s see which ASV has distance from root longer than 10.

long.branch <- dist.root[dist.root > 5]
long.branch
##    ASV601    ASV438 
## 20.819315  5.082379

Are these ASVs special?

dada2:::pfasta(refseq(ps)[names(long.branch)], ids = names(long.branch))
## >ASV601
## GGAGAGCCATGTCTTATAAGTATACCATACTCTTTTAACTTCTTAATGTTCTTTTCTATGCTTCGTATTCGAAAGAAATAAATTCACATTTTTTTCTGTGTATCAGTTAATTCATAATCCTGAAATCTCGGCTTCGATAGTCCATACTCTCTCTTATCCAGAAACACTGCACGATCTGTTCCCTTGAATACCGCACATTGTGTTTTAGCGAATGAAAAATAGTCCGATGTCAGCAACGCATAGGCATTCGTTGCCAGAAGCTGTCCCTCACTCTGCTTTAAAATCTTCCACCAACATC
## >ASV438
## CAAATCCGAACTGCAAAATACTCATGGAAAGAATCGCTCCTGTTTTTATTCCTTCTCCTAACATGGTAAGAACAGAAACTGGCACTCCAAATACAGAAAAACCACAAAGCATACAAATAAAAAGCCATCCGCCTCTCTGTGCCAGAACCTGGAAAAAATATTCTTTTCCAGAAACATCTTCCCCCGCAAAAGCACCCAAAAGGTAAAATACATGGACTGTTTCCTGTTTCCACTGCATCTTCCACAAAAGATTTGGAAGAACCGTTCCC

Result returned from BLAST (web interface):

ASV Taxonomy % identity
ASV601 Anaerobutyricum hallii 95.1
ASV438 uncultured Blautia sp. 97.3

Now we remove this outlier ASVs

all.taxa <- taxa_names(ps)
keep.taxa <- all.taxa[!(all.taxa %in% names(long.branch))]
ps.sub <- ps %>% prune_taxa(keep.taxa, .)

Convert ps.sub to proportion

ps.sub <- ps.sub %>% transform_sample_counts(function(x) x / sum(x))
sample_sums(ps.sub) %>% summary()
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       1       1       1       1       1       1

Test for Effect

plot_ord <- function(data, factor, method, distance){
    data.ord <- ordinate(data, method = method, distance = distance)
    p <- plot_ordination(data, data.ord, color = factor)
    p <- p + stat_ellipse(type = "t",geom = "polygon",alpha = 0)
    p <- p + ggtitle(str_c(factor,method,distance, sep = ' - '))
    print(p)
}
permanova <- function(data, formula, method, permutations=1e4, strata = NULL, core = max.core){
    dist.matrix <- phyloseq::distance(data, method=method)
    df <- data.frame(sample_data(data))
    model <- as.formula(paste0('dist.matrix~', formula))
    if (!is_null(strata)) {strata <- df[,strata]}
    result <- adonis2(model,
                      data = df,
                      permutations=permutations,
                      strata = strata,
                      parallel = core)
    return(result)
}
permdisp <- function(data, group, method, permutations=1e4, pairwise = FALSE, core = max.core){
    dist.matrix <- phyloseq::distance(data, method=method)
    df <- data.frame(sample_data(data))
    beta.disp <- betadisper(dist.matrix, group = df[,group])
    result <- permutest(beta.disp, permutations = permutations, pairwise = pairwise, type = 'centroid')
    return(result)
}

Household Effect

permanova(ps.sub, 'Household', 'bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##           Df SumOfSqs      R2      F    Pr(>F)    
## Household 48  11.7188 0.68925 2.2643 9.999e-05 ***
## Residual  49   5.2834 0.31075                     
## Total     97  17.0021 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Epileptic Effect

Ordination

plot_ord(ps.sub, 'Epileptic.or.Control','MDS','bray')

plot_ord(ps.sub, 'Epileptic.or.Control','MDS','wunifrac')

plot_ord(ps.sub, 'Epileptic.or.Control','NMDS','bray')
## Run 0 stress 0.2084864 
## Run 1 stress 0.2041351 
## ... New best solution
## ... Procrustes: rmse 0.05361631  max resid 0.4049279 
## Run 2 stress 0.2035787 
## ... New best solution
## ... Procrustes: rmse 0.03136479  max resid 0.2746453 
## Run 3 stress 0.2029836 
## ... New best solution
## ... Procrustes: rmse 0.01769698  max resid 0.07861787 
## Run 4 stress 0.208685 
## Run 5 stress 0.2030752 
## ... Procrustes: rmse 0.004255481  max resid 0.02643523 
## Run 6 stress 0.2088515 
## Run 7 stress 0.2064062 
## Run 8 stress 0.2162702 
## Run 9 stress 0.2029817 
## ... New best solution
## ... Procrustes: rmse 0.0006812619  max resid 0.004905256 
## ... Similar to previous best
## Run 10 stress 0.2040754 
## Run 11 stress 0.2100049 
## Run 12 stress 0.2053367 
## Run 13 stress 0.2113234 
## Run 14 stress 0.2106682 
## Run 15 stress 0.2312135 
## Run 16 stress 0.2039773 
## Run 17 stress 0.2041273 
## Run 18 stress 0.2122018 
## Run 19 stress 0.209445 
## Run 20 stress 0.2047825 
## *** Best solution repeated 1 times

plot_ord(ps.sub, 'Epileptic.or.Control','NMDS','wunifrac')
## Run 0 stress 0.1726183 
## Run 1 stress 0.1740023 
## Run 2 stress 0.1762893 
## Run 3 stress 0.1718306 
## ... New best solution
## ... Procrustes: rmse 0.0235026  max resid 0.173255 
## Run 4 stress 0.1872898 
## Run 5 stress 0.1706236 
## ... New best solution
## ... Procrustes: rmse 0.01010811  max resid 0.09705579 
## Run 6 stress 0.1718317 
## Run 7 stress 0.1861904 
## Run 8 stress 0.1899837 
## Run 9 stress 0.1706237 
## ... Procrustes: rmse 5.413324e-05  max resid 0.0003045123 
## ... Similar to previous best
## Run 10 stress 0.1726201 
## Run 11 stress 0.1696998 
## ... New best solution
## ... Procrustes: rmse 0.02436172  max resid 0.1715287 
## Run 12 stress 0.1947747 
## Run 13 stress 0.1822626 
## Run 14 stress 0.1876531 
## Run 15 stress 0.1802945 
## Run 16 stress 0.1864107 
## Run 17 stress 0.1865221 
## Run 18 stress 0.1856708 
## Run 19 stress 0.1891249 
## Run 20 stress 0.185231 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##      1: no. of iterations >= maxit
##     19: stress ratio > sratmax

PERMANOVA Test

permanova(ps.sub, 'Epileptic.or.Control', 'bray', strata = 'Household')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks:  strata 
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##                      Df SumOfSqs      R2     F Pr(>F)
## Epileptic.or.Control  1   0.1499 0.00882 0.854  0.111
## Residual             96  16.8522 0.99118             
## Total                97  17.0021 1.00000
permanova(ps.sub, 'Epileptic.or.Control', 'wunifrac', strata = 'Household')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks:  strata 
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##                      Df SumOfSqs      R2      F Pr(>F)  
## Epileptic.or.Control  1   0.1149 0.01615 1.5756 0.0189 *
## Residual             96   6.9986 0.98385                
## Total                97   7.1135 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PERMDISP

permdisp(ps.sub, 'Epileptic.or.Control', 'bray')
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
## 
## Response: Distances
##           Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.00251 0.0025142 0.1605  10000 0.6924
## Residuals 96 1.50379 0.0156645
permdisp(ps.sub, 'Epileptic.or.Control', 'wunifrac')
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
## 
## Response: Distances
##           Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.02183 0.0218264 2.2935  10000   0.13
## Residuals 96 0.91358 0.0095165

Breed Effect

a table of # of breed. frequency table

(data.frame(sample_data(ps.sub))$Breed.Group..1.) %>% table()
## .
##          Herder              NA Pointer Spaniel       Retriever     Scent Hound 
##              20               9              26              26               5 
##     Sight Hound        Sled Dog         Terrier 
##               3               6               3
plot_ord(ps.sub, 'Breed.Group..1.','MDS','bray')
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

plot_ord(ps.sub, 'Breed.Group..1.','MDS','wunifrac')
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

plot_ord(ps.sub, 'Breed.Group..1.','NMDS','bray')
## Run 0 stress 0.2084864 
## Run 1 stress 0.2040108 
## ... New best solution
## ... Procrustes: rmse 0.05685435  max resid 0.4005662 
## Run 2 stress 0.2089403 
## Run 3 stress 0.2093198 
## Run 4 stress 0.20545 
## Run 5 stress 0.2173889 
## Run 6 stress 0.219423 
## Run 7 stress 0.203331 
## ... New best solution
## ... Procrustes: rmse 0.03779231  max resid 0.3175873 
## Run 8 stress 0.2047923 
## Run 9 stress 0.2042308 
## Run 10 stress 0.2047717 
## Run 11 stress 0.2049309 
## Run 12 stress 0.2037882 
## ... Procrustes: rmse 0.03544949  max resid 0.3059876 
## Run 13 stress 0.2047891 
## Run 14 stress 0.2049428 
## Run 15 stress 0.2041995 
## Run 16 stress 0.2030192 
## ... New best solution
## ... Procrustes: rmse 0.01241197  max resid 0.0799837 
## Run 17 stress 0.2041269 
## Run 18 stress 0.2029852 
## ... New best solution
## ... Procrustes: rmse 0.005615208  max resid 0.03165945 
## Run 19 stress 0.211464 
## Run 20 stress 0.2040137 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##      5: no. of iterations >= maxit
##     15: stress ratio > sratmax
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

plot_ord(ps.sub, 'Breed.Group..1.','NMDS','wunifrac')
## Run 0 stress 0.1726183 
## Run 1 stress 0.184764 
## Run 2 stress 0.1697108 
## ... New best solution
## ... Procrustes: rmse 0.01847867  max resid 0.1443311 
## Run 3 stress 0.1747189 
## Run 4 stress 0.1724669 
## Run 5 stress 0.1726183 
## Run 6 stress 0.1724713 
## Run 7 stress 0.1872728 
## Run 8 stress 0.1706105 
## Run 9 stress 0.1921029 
## Run 10 stress 0.1791179 
## Run 11 stress 0.1820495 
## Run 12 stress 0.1853757 
## Run 13 stress 0.1800669 
## Run 14 stress 0.1921885 
## Run 15 stress 0.1849801 
## Run 16 stress 0.1773687 
## Run 17 stress 0.1884582 
## Run 18 stress 0.1819956 
## Run 19 stress 0.1792786 
## Run 20 stress 0.1912135 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     19: stress ratio > sratmax
##      1: scale factor of the gradient < sfgrmin
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

PERMANOVA

# same breed group v.s. diff breed group in household
# check if data are paired within household
if (!identical(sample_data(ps)$Household[seq(1,98,2)],
               sample_data(ps)$Household[seq(1,98,2)+1])) stop('data is not paired')

same.breed.group <- table(sample_data(ps)$Household, sample_data(ps)$Breed.Group..1.) %>% 
    apply(1, function(x) any(x == 2)) # if dog's grouping are same with in each house
same.breed.group;sum(same.breed.group)
##     1    10    11    12    13    14    15    16    17    18    19     2    20 
## FALSE  TRUE  TRUE  TRUE FALSE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE 
##    21    22    23    24    25    26    27    28    29     3    30    31    32 
##  TRUE FALSE  TRUE FALSE  TRUE  TRUE FALSE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE 
##    33    34    35    36    38    39     4    40    41    42    43    44    45 
##  TRUE  TRUE  TRUE  TRUE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE 
##    47    48    49     5    50    51     6     7     8     9 
##  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE
## [1] 39

Here we see 39 out of 49 household have dogs in same breed group

permanova(ps.sub, 'Household + Breed.Group..1.', 'bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##                 Df SumOfSqs      R2      F    Pr(>F)    
## Household       48  11.7188 0.68925 2.3509 9.999e-05 ***
## Breed.Group..1.  6   0.8179 0.04811 1.3127   0.06899 .  
## Residual        43   4.4655 0.26264                     
## Total           97  17.0021 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permanova(ps.sub, 'Household + Breed.Group..1.', 'wunifrac')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##                 Df SumOfSqs      R2      F    Pr(>F)    
## Household       48   4.8365 0.67991 2.2564 9.999e-05 ***
## Breed.Group..1.  6   0.3568 0.05015 1.3316    0.1275    
## Residual        43   1.9202 0.26993                     
## Total           97   7.1135 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PERMANOVA Test by Epileptic and Control Group

Here we

ps.control <- ps.sub %>% subset_samples(Epileptic.or.Control == 'Control')
ps.epileptic <- ps.sub %>% subset_samples(Epileptic.or.Control == 'Epileptic')
perm.breed.bray.ctl <- permanova(ps.control, 'Breed.Group..1.', 'bray')
perm.breed.wunif.ctl <- permanova(ps.control, 'Breed.Group..1.', 'wunifrac')
perm.breed.bray.ctl;perm.breed.wunif.ctl
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##                 Df SumOfSqs      R2      F Pr(>F)  
## Breed.Group..1.  7   1.6962 0.19479 1.4169 0.0319 *
## Residual        41   7.0114 0.80521                
## Total           48   8.7076 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##                 Df SumOfSqs      R2      F Pr(>F)  
## Breed.Group..1.  7   0.8423 0.21377 1.5925  0.032 *
## Residual        41   3.0977 0.78623                
## Total           48   3.9400 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
perm.breed.bray.epi <- permanova(ps.epileptic, 'Breed.Group..1.', 'bray')
perm.breed.wunif.epi <- permanova(ps.epileptic, 'Breed.Group..1.', 'wunifrac')
perm.breed.bray.epi;perm.breed.wunif.epi
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##                 Df SumOfSqs      R2      F Pr(>F)
## Breed.Group..1.  7   1.2777 0.15688 1.0898 0.2615
## Residual        41   6.8669 0.84312              
## Total           48   8.1446 1.00000
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##                 Df SumOfSqs      R2      F Pr(>F)
## Breed.Group..1.  7  0.55214 0.18052 1.2902 0.1023
## Residual        41  2.50647 0.81948              
## Total           48  3.05861 1.00000

Combined p value

\[X^2=-2 \sum_{i=1}^k \ln \left(p_i\right)\]

bray.p <- c(perm.breed.bray.ctl$`Pr(>F)`[1], perm.breed.bray.epi$`Pr(>F)`[1])
# Calculate the combined p-value with fisher's method
pchisq( -2 * sum(log(bray.p)), df = 4, lower.tail = FALSE)
## [1] 0.04826188
wunif.p <- c(perm.breed.wunif.ctl$`Pr(>F)`[1], perm.breed.wunif.epi$`Pr(>F)`[1])
# Calculate the combined p-value with fisher's method
pchisq( -2 * sum(log(wunif.p)), df = 4, lower.tail = FALSE)
## [1] 0.02200095

With using the weighted-unifrac distance, there’s at least one of the null hypothesis is rejected.

Drug Effect

Here we are only focusing on epileptic dogs.

plot_ord(ps.epileptic, 'Pheno.Y.N','MDS','bray')

plot_ord(ps.epileptic, 'Pheno.Y.N','MDS','wunifrac')

plot_ord(ps.epileptic, 'Pheno.Y.N','NMDS','bray')
## Run 0 stress 0.2030888 
## Run 1 stress 0.2049861 
## Run 2 stress 0.2059649 
## Run 3 stress 0.2101223 
## Run 4 stress 0.2165297 
## Run 5 stress 0.1990618 
## ... New best solution
## ... Procrustes: rmse 0.0459343  max resid 0.1625879 
## Run 6 stress 0.1999106 
## Run 7 stress 0.2132579 
## Run 8 stress 0.2146445 
## Run 9 stress 0.2087338 
## Run 10 stress 0.2049879 
## Run 11 stress 0.2020433 
## Run 12 stress 0.2239229 
## Run 13 stress 0.2080408 
## Run 14 stress 0.2194199 
## Run 15 stress 0.2106283 
## Run 16 stress 0.2055968 
## Run 17 stress 0.2177038 
## Run 18 stress 0.2281718 
## Run 19 stress 0.206055 
## Run 20 stress 0.2180574 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress ratio > sratmax

plot_ord(ps.epileptic, 'Pheno.Y.N','NMDS','wunifrac')
## Run 0 stress 0.1857 
## Run 1 stress 0.1942337 
## Run 2 stress 0.1854942 
## ... New best solution
## ... Procrustes: rmse 0.07777616  max resid 0.237673 
## Run 3 stress 0.1873453 
## Run 4 stress 0.1874356 
## Run 5 stress 0.19423 
## Run 6 stress 0.1906979 
## Run 7 stress 0.1927119 
## Run 8 stress 0.1814438 
## ... New best solution
## ... Procrustes: rmse 0.0508053  max resid 0.2199367 
## Run 9 stress 0.1831429 
## Run 10 stress 0.1918678 
## Run 11 stress 0.1932858 
## Run 12 stress 0.4028392 
## Run 13 stress 0.193139 
## Run 14 stress 0.1884053 
## Run 15 stress 0.1857624 
## Run 16 stress 0.1943164 
## Run 17 stress 0.1937669 
## Run 18 stress 0.2090482 
## Run 19 stress 0.1830395 
## Run 20 stress 0.1894807 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress ratio > sratmax

PERMANOVA

permanova(ps.epileptic, 'Pheno.Y.N','bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##           Df SumOfSqs      R2      F Pr(>F)
## Pheno.Y.N  1   0.2090 0.02566 1.2379 0.1963
## Residual  47   7.9356 0.97434              
## Total     48   8.1446 1.00000
permanova(ps.epileptic, 'Pheno.Y.N','wunifrac')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##           Df SumOfSqs      R2      F Pr(>F)
## Pheno.Y.N  1  0.06486 0.02121 1.0182 0.4158
## Residual  47  2.99375 0.97879              
## Total     48  3.05861 1.00000

Dispersion

permdisp(ps.epileptic, 'Pheno.Y.N','bray')
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq     F N.Perm Pr(>F)
## Groups     1 0.01951 0.019507 1.437  10000 0.2338
## Residuals 47 0.63801 0.013575
permdisp(ps.epileptic, 'Pheno.Y.N','wunifrac')
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
## 
## Response: Distances
##           Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)  
## Groups     1 0.03162 0.0316221 4.2019  10000 0.0467 *
## Residuals 47 0.35370 0.0075256                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# within epileptic dogs
table(sample_data(ps)$Seizure.Freedom..Y.N., sample_data(ps)$Epileptic.or.Control)
##      
##       Control Epileptic
##   NA       49         3
##   No        0        32
##   Yes       0        14
table(sample_data(ps)$Seizure.Control..Satisfactory.Unsatisfactory., sample_data(ps)$Epileptic.or.Control)
##     
##      Control Epileptic
##   NA      49         3
##   S        0        33
##   US       0        13

Sex Effect

Does health condition related to dog’s sex? Here we conduct a Chi-squared test.

# chisq test on epi male dog v.s. female.
# not within health condition
tb <- table(sample_data(ps)$Sex, sample_data(ps)$Epileptic.or.Control)
tb
##    
##     Control Epileptic
##   F      33        25
##   M      16        24
chisq.test(tb)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tb
## X-squared = 2.0698, df = 1, p-value = 0.1502

Ordination

plot_ord(ps, 'Sex','MDS','bray')

plot_ord(ps, 'Sex','MDS','wunifrac')

plot_ord(ps, 'Sex','NMDS','bray')
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.2261393 
## Run 1 stress 0.2450121 
## Run 2 stress 0.2605166 
## Run 3 stress 0.2507737 
## Run 4 stress 0.2445078 
## Run 5 stress 0.2364009 
## Run 6 stress 0.2282095 
## Run 7 stress 0.229298 
## Run 8 stress 0.2365472 
## Run 9 stress 0.2371477 
## Run 10 stress 0.2342887 
## Run 11 stress 0.2296419 
## Run 12 stress 0.2547756 
## Run 13 stress 0.2597614 
## Run 14 stress 0.2361598 
## Run 15 stress 0.2432431 
## Run 16 stress 0.2491273 
## Run 17 stress 0.235258 
## Run 18 stress 0.2485772 
## Run 19 stress 0.2446727 
## Run 20 stress 0.2526808 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress ratio > sratmax

plot_ord(ps, 'Sex','NMDS','wunifrac')
## Run 0 stress 0.1724074 
## Run 1 stress 0.1810284 
## Run 2 stress 0.1705573 
## ... New best solution
## ... Procrustes: rmse 0.0224693  max resid 0.1611144 
## Run 3 stress 0.1808744 
## Run 4 stress 0.1971364 
## Run 5 stress 0.189983 
## Run 6 stress 0.1817739 
## Run 7 stress 0.1836256 
## Run 8 stress 0.1919375 
## Run 9 stress 0.1854364 
## Run 10 stress 0.1723312 
## Run 11 stress 0.1754331 
## Run 12 stress 0.1843778 
## Run 13 stress 0.1705571 
## ... New best solution
## ... Procrustes: rmse 0.0004753476  max resid 0.003322762 
## ... Similar to previous best
## Run 14 stress 0.1842893 
## Run 15 stress 0.172415 
## Run 16 stress 0.1790648 
## Run 17 stress 0.1864631 
## Run 18 stress 0.1763068 
## Run 19 stress 0.1778902 
## Run 20 stress 0.1879778 
## *** Best solution repeated 1 times

PERMANOVA

permanova(ps.sub, 'Household + Sex', 'bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##           Df SumOfSqs      R2      F    Pr(>F)    
## Household 48  11.7188 0.68925 2.2419 9.999e-05 ***
## Sex        1   0.0561 0.00330 0.5156    0.9684    
## Residual  48   5.2272 0.30745                     
## Total     97  17.0021 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permanova(ps.sub, 'Household + Sex', 'wunifrac')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##           Df SumOfSqs      R2      F    Pr(>F)    
## Household 48   4.8365 0.67991 2.1305 9.999e-05 ***
## Sex        1   0.0068 0.00096 0.1443    0.9985    
## Residual  48   2.2701 0.31913                     
## Total     97   7.1135 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PERMDISP

permdisp(ps.sub, 'Sex', 'bray')
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.01527 0.015267 0.9824  10000 0.3235
## Residuals 96 1.49198 0.015541
permdisp(ps.sub, 'Sex', 'wunifrac')
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.00068 0.000684 0.0672  10000 0.8027
## Residuals 96 0.97753 0.010183

Age Effect

ggplot(sam) +
    geom_line(aes(x = as.numeric(Household), y = as.numeric(Age..months.), group = Household)) + 
    geom_point(aes(x = as.numeric(Household), y = as.numeric(Age..months.), colour = Epileptic.or.Control)) +
    xlab('Household') + ylab('Age')

PERMANOVA

permanova(ps.sub, 'Household + Age..months.', 'bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##              Df SumOfSqs      R2      F    Pr(>F)    
## Household    48  11.7188 0.68925 2.3061 9.999e-05 ***
## Age..months. 19   2.1074 0.12395 1.0477    0.3502    
## Residual     30   3.1760 0.18680                     
## Total        97  17.0021 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permanova(ps.sub, 'Household + Age..months.', 'wunifrac')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
##              Df SumOfSqs      R2      F    Pr(>F)    
## Household    48   4.8365 0.67991 2.2162 9.999e-05 ***
## Age..months. 19   0.9130 0.12835 1.0569    0.3676    
## Residual     30   1.3639 0.19174                     
## Total        97   7.1135 1.00000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1